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The paper presents an analytical theory and numerical simulations of the dipolar response of hydrated pro- 
teins. The effective dielectric constant of the solvated protein, representing the average dipole moment induced 
at the protein by a uniform external field, shows a remarkable variation among the proteins studied by nu- 
merical simulations. It changes from 0.5 for ubiquitin to 640 for cytochrome c. The former value implies a 
negative dipolar susceptibility of ubiquitin, that is a dia-electric dipolar response and negative dielectrophore- 
sis. It means that a protein carrying an average dipole of ~ 240 D is expected to repel from the region of 
a stronger electric field. This outcome is the result of a negative cross-correlation between the protein and 
water dipoles, compensating for the positive variance of the protein dipole in the overall dipolar susceptibility. 
This phenomenon can therefore be characterized as overscreening of protein's dipole by the hydration shell. 
In contrast to the neutral ubiquitin, charged proteins studied here show para-electric dipolar response and 
positive dielectrophoresis. The protein-water dipolar cross-correlations are long-ranged, extending approxi- 
mately 2 nm from the protein surface into the bulk. A similar correlation length of about 1 nm is seen for the 
electrostatic potential produced by the hydration water inside the protein. The analysis of numerical sim- 
ulations suggests that the polarization of the protein-water interface is strongly affected by the distribution 
of the protein surface charge. This component of the protein dipolar response gains in importance for high 
frequencies, above the protein Debye peak, when the response of the protein dipole becomes dynamically 
arrested. The interface response found in simulations suggests a possibility of a positive increment of the 
high-frequency dielectric constant of the solution compared to the dielectric constant of the solvent. This 
analysis provides a theoretical foundation for experimentally observed positive increments of the absorption 
of THz radiation by protein solutions. 

Keywords: Protein solvation, dielectric response, dielectrophoresis, protein electrostatics, THz absorption, 
cavity field 



I. INTRODUCTION 

Polarization of the interface is an important compo- 
nent of the response of a polar substance to an external 
field. The standard approach of Maxwell's electrostat- 
ics assumes that the interface cuts through the polarized 
dipoles of the dielectric, leaving their monopoles at the 
surface (Fig. [1^) • The density of these monopoles is the 
surface charge density ap. It is given by the projection 
Pn of the dipolar polarization vector P on the outward 
normal n to a continuous dielectric mediumfii^ This sur- 
face charge is opposite in sign to an external charge and 
so the field of the surface charges compensates the exter- 
nal field. The sum of the two fields makes the Maxwell 
field inside the dielectric, which is lower in intensity than 
the external field. 

The same basic considerations apply to the problem of 
solutions polarized by a uniform external field Eq. The 
interface is now a closed surface enveloping each solute. 
The polarization field, uniform in the bulk, becomes in- 
homogeneous close to the solute-solvent interface. It gen- 
erates positive and negative lobes of the surface charge 
density (Fig. [2|) integrating into an overall interface (sub- 
script "int") dipole Mq"'. Its calculation is generally a 
complex problem involving both the effects of the so- 



lute shape and the alteration of the liquid structure by 
the solute surface multipoles. A closed-form solution is, 
however, possible in the framework of standard dielec- 
tric theories for a spherical void in a dielectric^ when all 
specifics of the solute-solvent interactions are neglected 



Mj,"* = -3r!oP/(2e, 



(1) 
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Here, Eg is the solvent dielectric constant, Hq is the solute 
volume, and subscript "0" is assigned throughout below 
to the solute parameters. The orientation of the inter- 
face dipole is opposite to the uniform polarization of the 
medium P. 

The appearance of Mq"' is an interfacial, but not nec- 
essarily a surface phenomenon. This interface dipole is 
the integral effect of the inhomogencous polarization sur- 
rounding an excluded volume of the solute. This polar- 
ization perturbation in fact propagates quite far into the 
bulk, as we show below, and can be taken fully into ac- 
count only in the thermodynamic limit for the solution, 
which we represent below as the fc — > limit in the in- 
verted Fourier space of wavcvectors k. The surface charge 
density ap is just a convenient mathematical represen- 
tation of this highly non-local physical reality in terms 
of properties assigned to an infinitely thin mathematical 
surface enveloping the solute. 

The interface dipole arising from a void in a uniformly 
polarized liquid will in turn polarize the surrounding sol- 
vent. As a result, the dipole moment of a uniformly po- 
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FIG. 1. Cartoon of two extreme polarization patterns of a 
polar solvent at the surface of a spherical solute with a positive 
charge at the center. In panel (a), the solvent dipoles align 
with the electric field of the solute charge. This alignment 
results in surface charge density ap = Pn < 0. In panel 
(b), the surface dipoles preserve their preferential in- plane 
orientations characteristic of a free surface of a polar liquid.'' 
This orientational pattern produces no surface charge, ap = 
0. 

larized solution is given as 

M, = M"^ - rioP - (2/3) (e,, - l)Mj,"*. (2) 

The first summand in this equation is the dipole moment 
of a uniformly polarized homogeneous liquid. The seeond 
term is the dipole moment reduced from the liquid by 
putting a solute of volume JIq into it. Finally, the last 
summand is the polarization of the liquid by the interface 
dipole. 

The interface dipole M™* exists at any closed interface 
in a polarized medium, even in the absence of solute's 
own charges. When calculated according to Eqs. ([T|) and 
([2|). it will lower the dielectric constant of the solution 
e compared to the dielectric constant of the homoge- 
neous liquid. In contrast, if the solute possesses its own 
dipole, it will align along the external field Eq producing 
an average permanent dipole (Mo)_b. This dipole mo- 
ment will enhance the dielectric response of the solution. 
The overall dipole moment associated with a solute will 
be the sum of the intrinsic permanent dipole and the 
dipole induced at the dielectric interface^ 

Mo = (Mo)ij + Mi,'^*, (3) 

where (. . .)e refers to the statistical average in the pres- 
ence of the field. 

Both the permanent and interface components of Mg 
depend on the ability of the interface to polarize. The 
very basic physics outlined in Fig.[T^ assumes the dipoles 
of the medium to have the ability to freely change their 
orientations in order to align along an external electric 
field. While this is probably the case for solvation of small 
ions in polar liquids, solvation of larger solutes might 
present a challenge to this picture. 

Dipoles of polar liquids preferentially orient in-plane at 
interfaces4r— Unless an external field rotates the dipoles 
off-plane, such orientational structure eliminates the sur- 
face charge since crp = P„ — in this case (Fig. [Us). The 




FIG. 2. Cartoon of the interface polarization of a spheri- 
cal void in a uniformly polarized polar liquid. The surface 
charges at the interface produce negative and positive lobes 
of the overall surface charge density responsible for the inter- 
face dipole Mo"* = xi^^o^-o- When the surface charge density 
disappears because of the in-plane alignment of the surface 
dipoles, the dipole subtracted from the solution is simply the 
product of the solute volume and the uniform polarization of 
the bulk [Eq. ((2|]. This is the Lorentz scenario of the inter- 
face polarization corresponding to xi = 0. In contrast, the 
Maxwell scenario anticipates a non-zero surface charge den- 
sity, which reduces the dipole taken from the solution from 
the Lorentz value — P^fJo to the Maxwell value given by Eq. 
((!]). This scenario also anticipates a nonzero and negative 
dipolar interface susceptibility xi given by Eq. (fTS)) . 

standard boundary conditions of continuum electrostat- 
ics do not apply to this scenario which implies 

Mj,"' =0, (4) 

instead of the standard electrostatic result listed in Eq. 

Equations ([2]) and ([4]) suggest that carving a void from 
a dielectric removes the dipole moment equal to the prod- 
uct of the uniform bulk polarization and the void volume. 
Such a solution would appear in the standard theories of 
dielectrics^^ if the void in the dielectric had the inter- 
face of a Lorentz virtual cavity, which has no surface 
polarization by definition. We will dubb this outcome, 
corresponding to ap = in Fig.[T|D, as the "Lorentz sce- 
nario". On the other hand, when Mq"' from Eq. ^ is 
substituted into Eq. ([2]), the two last terms combine into 
the same M"*, which becomes the dipole subtracted from 
the homogeneous liquid upon insertion of the solute. The 
rules of electrostatics therefore predict that the dipole re- 
moved from the solution will be, at 3> 1, just a small 
fraction of the Lorentz dipole. Since this scenario follows 
from the standard electrostatics with Maxwell's bound- 
ary conditions at the dielectric interface, we will call this 
outcome, and the corresponding interface polarization, 
the "Maxwell scenario" (Figs. [TJa, and [2|). 

Water presents a particularly important study case for 
the Lorentz scenario of Eq. Large solutes, over 1 

nm in size,— break the network of hydrogen bonds of 
bulk water, resulting in preferential in-plane orientation 
of the interfacial water dipoles4i^— This phenomenon, 
general for polar fiuids^iii^ is further amplified for water 
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interfaces by the high energy of water's hydrogen bonds<S 
The orientational structure of interfacial water has 
macroscopicaUy observable consequences, both mechan- 
ical and electrostatic. Mechanical consequences include 
the rotation of hydrated nanometer solutes by external 
fieldsi^ and slipping of the boundary layers in the hydro- 
dynamic flow.— For electrostatic observables, the electric 
field inside cavities formed in uniformly polarized dipolar 
liquids^ and in water— seem to follow the scenario of an 
unpolarized interface sketches in Fig. [TJj. One can there- 
fore anticipate that in-plane dipolar orientations might 
be preserved, at least in patches, if the external field is 
not sufficiently strong to compete with interface hydro- 
gen bonds forcing water dipoles orient in-plane. If this 
is the case, the boundary conditions of the dielectric re- 
sponse problem will alter, thus affecting all relevant polar 
response functions. 

Electrostatic interactions are critical for biological 
functionJ^"— Most biomolecules and all cellular mem- 
branes carry charges Electrostatic solvation and inter- 
actions affect the stability of folded proteins and their 
aggregation and crystallization^i"— Therefore, the ques- 
tion of what is the dipolar polarization at the interface of 
a hydrated biomolecule is critical for structural biology 
and bioenergeticsi^ 

Surfaces of proteins, and of all biomaterials in 
general, arc obviously chemically and electrostatically 
heterogeneous 1^ One therefore cannot expect a clear-cut 
scenario of either in-pane dipoles or dipoles fully aligned 
along the electric field. This study in fact shows that 
none of the scenarios sketched in Fig. [1] presents a com- 
plete description of hydrated proteins, which explore a 
much wider range of possibilities allowing them to tune 
their response to global and local in vitro fields. To grasp 
this complexity, we ask what would be a minimal set of 
coarse-grained parameters describing the water-protein 
interface. We approach this question here by first pre- 
senting an analytical theory framing the problem in terms 
of a set of interface susceptibilities, followed by numerical 
simulations of several hydrated globular proteins. The 
property of interest is the average dipole moment Mq 
[Eq. induced at the protein by an external electric 
field. As such, this is a fundamental and well-defined 
physical quantity related to broad-band dielectric spec- 
troscopy of solutions, not considered here, and directly 
probed by dielectrophoresis of protein solutions^^"— dis- 
cussed below. 



II. DIELECTRIC CONSIDERATIONS 

According to the separation of the solute dipole into 
the intrinsic permanent and interfacial components, one 
can define linear dipolar susceptibilities for the corre- 
sponding dipoles along the external field. If the 2-axis of 
the laboratory frame is set along the external field, one 
gets for the dipolar ( "d" ) and interfacial ( "int" ) parts of 



the response 

x^' = M-^:/{n,E,). 

We first focus on the solute permanent dipole and con- 
sider the first-order perturbation expansion for the cor- 
responding susceptibility^ 

xg = (/3/31)o)('5Mo-5M), (6) 

where j3 = l/{k^T) is the inverse temperature. Further, 
(5Mo = Mo -(Mo) and 5M. = M-(M) arc the deviations 
of the solute (Mo) and total sample (M) dipole moments 
from their corresponding average values; (...) refers to 
a statistical average in the absence of an external field. 
Since the solute can in principle be charged, keeping the 
variation (JMq eliminates the dependence of its dipole 
on the origin of the system of coordinates. On the other 
hand, (5M can be replaced with M if the sample is neutral 
and isotropic. 

The dipolar susceptibility 

Xo = x[? + xi)"* (7) 

connects average dipole of the solute to a weak external 
field Eq which varies on a length-scale significantly larger 
than the solute dimension. This electric field creates an 
excess chemical potential of the solute^ 

A/io = -\mo, Eo, (8) 

where Mq is given by Eq. ([3]). The dipolar susceptibil- 
ity therefore follows from the derivative of A^o over the 
external field strength: xqUqEq — —dAfiQ/dEQ. Since 
the external field Eq is well defined by the charge density 
at the plates of a planar capacitor in the dielectric ex- 
periment or by the light intensity in absorption measure- 
ments, the corresponding susceptibility is a well defined 
parameter as well. 

This susceptibility is obviously distinct from the sus- 
ceptibility of the protein solution responding to the 
macroscopic Maxwell field E. One can consider the av- 
erage solute dipole created in response to E. In that 
case, one needs a connection between the two fields. This 
connection, E = E^/e (where e is the solution dielectric 
constant), is again straightforward in dielectric measure- 
ments, but depends on the simulation protocol in numer- 
ical simulations 1^1^ It simplifies, however, significantly 
for tin-foil implementation of the Ewald sums represent- 
ing Coulomb interactions in simulations with periodically 
replicated simulation celli^ In that case, which is fol- 
lowed in this study, E = Eq and Xo from MD trajecto- 
ries gives the response to the Maxwell field E. A more 
elaborate theory, which we present below, is, however, 
needed to obtain the interface susceptibility x"*. Once 
this is done, one can follow the standard convention to 
introduce the statistical (that is obtained from the vari- 
ance of the dipole moment) dielectric constant^^ of the 
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protein 

eo = l+47rxo- (9) 

The dielectric constant eo , used here to quantify the so- 
lute dipolar response, is neither the dielectric constant of 
the solution e reported by the dielectric experiment^"— 
nor it is the screening parameter used to describe 
Coulomb interactions between charges inside the solute, 
such as atomic charges of protein residues j^^^i^ We also 
stress that eg defined by Eqs. ([7|) and ^ represents the 
dipolar response of the entire protein (irrespective of its 
shape) and tells nothing about dielectric properties of 
any region inside the protciui'^S The definition of the so- 
lute dipolar response in terms of the external field Eq, 
instead of the Maxwell field E, is convenient for a num- 
ber of reasons, including the fact that one does not need 
to calculate the solution dielectric constant e, which re- 
quires additional modeling. 

The susceptibility 

Xo = Xoo + Xos (10) 

is a sum of the direct solute component xoo cx (SMq) 
and a cross-correlation term xos ((5Mo • (5Ms), where 
Ms = M — Mo is the solvent dipole moment. While 
Xoo I the dielectric susceptibility defined for a finite solute 
volume ilor^ is obviously positive, one wonders what is 
the sign and relative magnitude of the xos component)^ 
The answer depends strongly on the details of the dipolar 
polarization of the solute-solvent interface. 

One needs to recognize that the correlator of the solute 
dipole Mo with the entire sample dipole M in Eq. 
instead of the self-correlator (SMq) typically appearing 
in macroscopic theories of dielectrics, substitutes for the 
boundary conditions used in these theories. The bound- 
ary conditions represent the physical fact that any treat- 
ment of the dipolar polarization of a finite sample should 
include surface charges if the sample is placed in vacuum 
or the polarization of the surrounding medium if the fi- 
nite sample is a part of an infinite dielectric materiali^ 
This notion also implies that understanding and po- 
tentially modeling of the cross-correlation susceptibility 
Xos oc ((5Mo • SlsAg) allows one to substitute the bound- 
ary conditions of standard electrostatics, which rely on 
the bulk dielectric constant jii^ with microscopic rules in- 
corporating various polarization scenarios, such as two 
extremes sketched in Fig. [TJ This perspective is particu- 
larly important for studies of solvation at the nano-meter 
scale. The number of first-shell waters of a typical globu- 
lar protein reaches the magnitude of A^i ~ 200 — 500, sep- 
arating them, and potentially other nearby shells, into a 
sub-ensemble. The properties of this sub-ensemble might 
dramatically differ from those of the bulk solvent. The 
language of intcrfacial susceptibilities might better grasp 
this reality than bulk properties typically used in stan- 
dard theories to construct the interfacial response func- 
tions. 

Accordingly, we introduce below the dipolar suscepti- 
bility of the surface charge density xi which will incor- 



porate all possible boundary conditions at the surface of 
the solute and will reproduce the Lorentz and Maxwell 
scenarios as special cases. We will use this susceptibility 
to both provide a connection between xos and xoo com- 
ponents of the dipolar response function Xo and derive a 
relation for Xo"' ■ This formalism will allow us to analyze 
the results of numerical simulations of protein solutions 
presented next. 

III. RESPONSE FUNCTIONS 

A general insight into the understanding of the sol- 
vent polarization surrounding a solute can be gained from 
the inhomogcneous response function of the dipolar po- 
larization in the solute's vicinity. These types of prob- 
lems typically involve integral convolutions in real space, 
which become algebraic relations in inverted k-spacc. In 
the presence of the solute, the dipolar response function 
x(ki,k2), which is a rank-two tensor, becomes a func- 
tion of two wave- vectors, ki and k2.— This is a reflection 
of the inhomogcneous character of the problem, in con- 
trast to the response function of the homogeneous solvent 
Xs(k), which depends on one wave- vector only. 

The modification introduced to the solvent response by 
the solute can be given as an inhomogcneous correction 
to Xs(k) as follows 

X(ki, k2) = Xs(ki)5ki,k. - (xi/xf )Xo(ki, k2). (11) 

Here, the Kronecker delta ^ki.k2 is normalized to the sam- 
ple volume V, Sq q = V. Further, the homogeneous liq- 
uid has isotropic symmetry, which is broken in k-space 
by the wave-vector k — k./k introducing axial symmetry 
to the problem. The second-rank tensor functions are 
then fully described by two scalar projections, longitu- 
dinal (L) and transverse (T), onto two diadic tensors,— 
= kk and 3^ = 1 — kk. Correspondingly, the re- 
sponse function of the homogeneous solvent is given as 
the sum of the longitudinal and transverse components 
;^^(k) = XsW^^+X^{k)3'^. Their fc = values are bulk 
susceptibilities connected to the liquid dielectric constant 

Xf (0) =(e. - l)/(47re,), 
Xj(0)=(6.-l)/(47r). 

The parameter xi in fi'ont of heterogeneous part of 
the response xo(ki,k2) in Eq. (jlip is the dipolar sus- 
ceptibility of the solute interface coarse-grained into a 
spherical surface. It appears from the expansion of the 
surface charge density induced by a uniform external field 
in Legcndre polynomials of the polar angle 6 between 
a radius-vector at the surface and the external field, 
ap{6) = J2e'^ePe{cos9). Only the first-order, dipolar 
component £ = 1 contributes to the dipolar polarization 
field P(r) of the liquid. The susceptibility xi connects 
(Ti to the field strength 

0-1 = XiEq. (13) 
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Correspondingly, the interface dipole is calculated by 
multiplying <7p{9) with the surface radius- vector and in- 
tegrating over the closed surface. The result is 

APof = Xi^io^^o (14) 

and, from Eq. Xq"' = xi- 

The derivation of x(ki,k2) for a void in a polar liq- 
uid that we discuss in Appendix |^ yields the standard 
Maxwell result for Af™* in Eq. ([1]) and a negative xi 

Xf = Xf (0). (15) 

In contrast, the Lorentz scenario corresponds to xi = 0. 
Since we do not want to be limited by the solution de- 
scribing a void and instead want to introduce an inter- 
facial polarization induced by the solute, xi in Eq- (fTTj) 
is an unspecified susceptibility calculated below from nu- 
merical simulations. 

The polarization of the solvent in response to a field 
Eo(r), which might include both solute and external 
charges, is given by the convolution in inverted space 

P(k) =x(k,k')*Eo(k'), (16) 

where tildes over vectors specify inverted-space fields and 
asterisks denotes both the tensor contraction and k-space 
integration, such as 

A*B = ■ B dk/(27r)^ (17) 

We now want to apply the inhomogeneous response 
function in Eq. (fTTj) to calculate the cross-correlation 
term xos in Eqs. © and (|10p . In order to approach 
this problem, we will calculate the solvent dipole in- 
duced by the inhomogeneous field of an instantaneous 
solute dipole Mq. We will assume a certain separation 
of time-scales to perform this calculation. Specifically, 
the solvent is assumed to be much faster than the solute 
and thus able to follow adiabatically every instantaneous 
configuration of the solute electric field. This approxima- 
tion is typically correct for hydrated proteins since high- 
frequency protein vibrations produce a relatively minor 
effect on the protein field sensed by hydration water 

For this problem, the external electric field now be- 
comes 

Eo(k) = T(k)-Mo, (18) 

where 

T(k) = -4.^D, (19) 

is the k-space dipolar tensor of a spherical solute with an 
effective radius i?; D = 3kk — 1 = 2J^ — and jtix) is 
the spherical Bessel function. 

The overall solvent dipole can be obtained by in- 
tegrating the dipolar polarization field P(r) over the vol- 
ume occupied by the solvent, ~ V — flo- This direct- 
space integration is equivalent, in inverted space, to tak- 
ing fc = limit for the polarization of the entire sample 



solvent solvent 



(a) I (b) I 
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FIG. 3. Cartoon of the solvent polarization around a dipolar 
solute. A higher number of solvent dipoles in the equatorial 
plane of the solute dipole Mq (a) yields a negative value of 
the response function xoa. In panel (b), the addition of equa- 
torial charges to the solute does not alter the solute dipole 
calculated relative to the center of the sphere, but polarizes 
equatorial waters in a direction orthogonal to Mq. The re- 
sults is a positive xos , as indeed found for all charged proteins 
here. 

(1^ — > oo in the thermodynamic limit) and subtracting 
the convolution of the polarization field with the step 
function 0o(r), equal to unity within the solute and zero 
everywhere elsci^ One then gets for the inverted-space 
fields 

M, = P(0) - P(k) * 9o{k), (20) 
where P(k) is given by Eq. (fT6|) and 

00 (k) = / e^'^'-dr. (21) 

We show in Appendix El that the second summand 
in Eq. pO]) is zero when the continuum, fc — >■ limit 
is taken in the homogeneous solvent response function 
Xs(k). We will use this limit throughout below since 
the typical size of the protein 2R significantly exceeds 
the diameter of water. Within this approximation one 
arrives (see Appendix |X| at xos and then at the following 
connection between the xos a-nd xoo response functions 
in terms of the interfacial dipolar susceptibility Xi 

Xos = -^^^(l-(4V3)e.Xi)Xoo. (22) 

3 Es 

From this equation and Eq. (|10p . one additionally have 

2^ = iL±l + ?I,„-lto. (23) 

Xoo 3es 9 

If the term in the brackets in Eq. (|22p is positive, as 
is the case in both the Lorentz and Maxwell scenarios, 
Xos is negative. Figure [3^ illustrates the physical ori- 
gin of this result. The solvent polarization is a sum of 
two compensating contributions. The solvent dipoles at 
the poles of the solute dipole will predominantly orient 
parallel to Mq, while equatorial solvent dipoles will ori- 
ent antiparallel to Mq. Since there are more equatorial 
dipoles than there are pole dipoles, xos < for the overall 
solvent polarization. 
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It is instructive to see what are the numerical outputs 
for Xo in the MaxweU and Lorentz scenarios sketched in 
Fig. [2 In the Lorentz case of xi = one gets 



Xo 



+ 2 

3es 



-Xoo- 



(24) 



The correction term in front of xoo is the Lorentz cav- 
ity field which indeed appears in polar response when 
the surface of a cavity cut from the dielectric is not 
polarized i^ii^ 

When the Maxwell result [Eq. (flSl) ] for the susceptibil- 
ity xi is used in Eq. ([^ one gets 



2e, 



1 



Xoo- 



(25) 



This correction factor is the well-known cavity field of the 
theories of dielectrics i^ii^ These arguments stress again 
(see above) that specifying an algorithm of calculating 
Ms in terms of Mq leads to a route to formulate a theory 
of polar response, with the prescription of Figs.[T^ and [2] 
corresponding to Maxwell's electrostatics i^'^'^° 

Equations (|24|) and (|25|) are special cases of a more 
general result 



Xo 
Xoo 



(26) 



connecting the cavity field inside the solute Ec with the 
ratio of two response functions. The notion "cavity" here 
implies that this field is produced by the solvent polarized 
by the external field and does not include any reaction 
fields of the solute charges. While the two contributions 
into the overall electric field inside the solute might seem 
to be hopelessly entangled, they are in fact separable in 
the frequency domain. The solute and corresponding re- 
action fields are dynamically frozen at frequencies above 
the Debye peak of the solute, in most practical cases well 
below the Debye peak of the water component of the so- 
lution. The observation of the dielectric increment of the 
water's Debye peak as a function of the solute concentra- 
tion in the limit of ideal solution provides a direct access 
to the cavity field^ and, by Eqs. ([23]) and (l26|), to the 
ratio of the two susceptibilities and xi- The relation of 
these susceptibilities to dielectrophoresis discussed below 
is yet another connection to the laboratory experiment. 

A significant difference in the solute dipolar response 
predicted by Lorentz and Maxwell polarization scenar- 
ios offers an opportunity to control the interaction of the 
solute dipolc with an external field by changing the dis- 
tribution of the surface charge. For example, consider 
the effect of surface charges placed in the equatorial re- 
gion of the global solute dipolc in Fig. [3)d. These charges 
will not change the overall solute dipole calculated rela- 
tive to the sphere's center, but will orient solvent dipoles 
perpendicular to the global dipolar field of the solute and 
can potentially invert the sign of xos • This is indeed what 
we find in our simulations of charged hydratcd proteins. 



One can also anticipate some combination of the solute 
shape and charge distribution that will produce a nega- 
tive net result for xoj when negative xos exceeds in mag- 
nitude Xoo- This outcome, which can be characterized 
as overscreening^ of the solute dipole by the hydration 
layer, would correspond to a dia-electric response of the 
solute, i.e. its repulsion from the region of a more intense 
electric field. We find this result in our simulations of the 
neutral ubiquitin (ubiq) protein. 



IV. RESULTS OF MD SIMULATIONS 



Equation (|22|) is a central result of our derivation. It 
connects the two component of the solute dipolar sus- 
ceptibility and the interface moment A/q"*^ to one single 
property, the dipolar susceptibility of the interface xi- 
This susceptibility is a coarse-grained parameter incor- 
porating the effects of both the shape and surface charge 
distribution of the solute. Since both xos and xoo SuTC in 
principle accessible from numerical simulations, Eq. (j22p 
allows us to construct the overall dipolar response of a 
hydrated solute without additional assumptions on the 
nature of polarization boundary conditions at the inter- 
face. 

Four globular proteins, reduced form of cytochrome 
c (cytC, 100 ns, PBD database 2B4Z), ubiquitin (ubiq, 
172 ns, lUBQ), lysozyme (lys, 153 ns, 3FE0), and re- 
duced form of cytochrome B562 (cytB, 123 ns, 256B) 
were studied by long, 100-172 ns all-atom Molecular Dy- 
namics (MD) simulations. The overall length of the sim- 
ulation trajectories was > 0.5 /iS. All proteins were sol- 
vated in large numbers of TIP3P waters to achieve mM 
protein concentration typically used in experimental so- 
lution measurements. The number of waters in the simu- 
lation box were N, = 33189 (cytC), 27918 (ubiq), 27673 
(lys), and 33268 (cytB). From these proteins, cytB is the 
only protein with the overall negative charge Z (Table [l|, 
and it was chosen mostly for that reason. 

The importance of long simulations has been rec- 
ognized in previous attempts to address the dielectric 
properties of protein solutions^^— but they had never 
reached the length and system size reported here. In par- 
ticular, Steinhauser and co-workers^^i^i^ have pointed 
to a significant contribution of Xos susceptibility to the 
overall dipolar response of the solution, but could reach 
only qualitative conclusions from trajectories available to 
them (~ 13 — 30 ns). Many previous attempts to sim- 
ulate protein solutions had suffered from even shorter 
trajectories and far smaller simulation boxes. The aim 
of this round of simulations is to extend previous simula- 
tion studies to reliably calculate both xoo a-nd xos com- 
ponents of the solute dipolar response. We therefore do 
not approach here the calculation of the overall dielectric 
response of the protein solution and instead focus on the 
question of how polar, as probed by xoj a hydratcd pro- 
tein can be. The details of the simulation protocol can 
be found in the Supplementary Material (SM)i^ and we 
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proceed directly to the results. 



A. Static properties 

The proteins studied here are similar in their effective 
size, as calculated from the volume inside the solvent- 
accessible surface enveloping the molecule.— The average 
dipole moments (Mq) calculated from the protein par- 
tial charges and averaged over the trajectories are also 
close in magnitude and are in basic agreement with dipole 
moments typically reported by dielectric spectroscopy of 
protein solutions^ft (Table |T| . The values of the dipole 
moments were calculated here relative to the protein cen- 
ter of mass since these dipoles appear in the equations 
of motion establishing the torque applied to the protein 
by an external electric field4i While the protein size and 
the dipole moment appear to be generic for the set of 
globular proteins studied here, the dipolar susceptibility, 
i.e. the variance of the protein dipole is quite specific for 
a given protein. We find a remarkably broad range of 
dipolar susceptibilities xo and the corresponding values 
of eo = 1 + inxo among the proteins studied here (Table 

The polarity of the protein, as measured by xo or eo, 
does not seem to correlate with either the magnitude or 
the sign of the overall protein charge. In fact, values of 
eg are close for positively charged lys {Z = 4-7) and neg- 
atively charged cytB {Z — —5). A significant variation 
in the values of eq found here (Table |T| seems to orig- 
inate from differences in surface charge distributions of 
the proteins and the corresponding polarization of the 
hydration water. 

The cross-correlation susceptibility xos also varies sig- 
nificantly among the proteins, both in magnitude and 
sign. We find xos positive for the charged proteins. This 
observation can be explained in terms of the water po- 
larization scenario pictured in Fig. [SJd. In contrast, xos 
is negative for the neutral ubiq, in agreement with the 
dielectric arguments presented above. However, in con- 
trast to the standard expectations, Xo is negative. This 
outcome implies a dia-electric response or negative di- 
electrophoresis. The relative error of this claim is 20% 
(Table H]). This is because Xo comes as a small num- 
ber from subtraction of two relatively large numbers, xoo 
and xos, and converges very slowly even on the longest 
trajectory we have run in this study (see the SM^). 

The results for the interface susceptibility xi are con- 
sistent with the results for xos (Table |l|. The suscepti- 
bility xi not only exceeds the Maxwell Xi^ [Eq. in 
magnitude, but has the sign opposite to it. A positive xi 
physically means that the dipole produced by the surface 
charge density ap (Fig. [2]) is oriented along the field and 
not opposite to the field as the Maxwell scenario would 
suggest. As for xos, the origin of this outcome should 
be sought along the lines illustrated in Fig. which 
shows that different regions of the interface polarization 
can add constructively or destructively depending on the 




5 10 15 20 25 30 



r/A 

FIG. 4. Reduced functions Xos{r)/xos (circles) and N{r)/Ns 
(squares) calculated from the layer of water of thickness r 
surrounding the vdW surface of cytC. Here, xos is calculated 
for the entire simulation cell containing Ns waters. The solid 
lines is the fit of Xos{r)/xos to an exponential function, 1 — 
exp[-(r - a)/A], with A = 24 A and a = 2.6 A. The dashed 
line connects the simulation points. 

distribution of the surface charge. 

Spatial correlation between the protein and water 
dipoles are found to be long-ranged. This is illustrated 
in Fig.[3]wherc wc plot Xos(^) and N(r) reduced to their 
values obtained from the entire simulation box. Function 
N{r) is a number of waters in a shell of thickness r as 
measured from the protein van der Waals (vdW) surface; 
Xos{r) oc ((5Mo ■5Ms(r)) is the cross-correlation function 
obtained from the total dipole moment Ms(r) of all wa- 
ters within the same shell. Although Xos('') clearly goes 
faster to its saturation limit than the number of waters 
in the cell, it reaches only half of its cell value for seven 
water shells around the protein. In fact when xos (r) /xos 
is fitted to an exponential decay function (solid line in 
Fig. |4]), the corresponding correlation length turns out 
to be 24 A. 

The nearly disappearance of the cross-correlation be- 
tween the protein and shell dipoles at the shell thickness 
r approaching the limit of one hydration layer implies 
that waters in the first solvation shell do not much corre- 
late with the overall protein dipole and are more driven 
by the local fields and vdW forces of the surface groups. 
It is only more distant layers that correlate more exten- 
sively with the global electrostatics of the solution. This 
picture contrasts with a more localized density response 
of the interface shown in Fig. [5l The density profile 
p{r) = N{r)/n{r) {ft{r) is the shell volume) peaks at the 
interface, indicating wetting of the protein surface by wa- 
ter, and then decays to approximately the bulk density 
within 3 solvation layers (~ 10 A). 

We also show in Fig. [5] the dielectric susceptibility of 
the surface waters determined by correlating the dipole 
moment of the shell Ms(r) with the total dipole moment 
of water Mg.— We define two susceptibilities, normalized 
to the number of shell waters, XNii^) oc 7V(r)~^ ((5Ms(r) ■ 
(5Ms), and the susceptibility normalized to the shell vol- 
ume, xn{r) oc n{T)~^ {S'M.sir) ■ S'M.g)- While the former 
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TABLE I. Dielectric parameters of hydrated proteins from MD simulations. 



Protein^ 


Z 


i?(A) 


eo 


(Mo) (D) 


d b 
Xo _ 


Xoo 




xi."' X 10^ 


xi/xf 




Lys 


+7 


19.2 


62 


223 


4.9 


3.9 


1.0 


3.7 


-2.8 


0.3 


CytC 


+6 


18.7 


639 


367 


51 


38 


13 


3.9 


-3.1 


5.7 


Ubiq 





17.8 


0.5 


244 


-0.04 


0.17 


-0.21 


-2.3 


1.78 


0.2 


CytB 


-5 


23.5 


96 


196 


7.5 


6.1 


1.4 


3.6 


-2.8 


21.6 



^ Z is the overall charge of the protein, R is the effective radius, the dipole moment (Mo) is calculated from the protein charges relative 
to the center of mass and averaged over the simulation trajectory. 

The standard deviation of Xq calculated according to the standard procedures explained in more detail in the SMfi are: 0.1(2%) (lys), 
1.6(3%) (cytC), 6 X 10-3(20%) (ubiq), 0.3(2%) (cytB). 

Compressibility of the first hydration shell, fti = {{SNi)^) / (Ni) , Ni is the number of waters in the first shell defined as the layer of 3 
A thickness measured from the protein vdW surface. 
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FIG. 5. Dependence on the thickness r of the solvation shell 
of the water number density p(r) = N{r)/Q{r) and two di- 
electric susceptibilities, XfiC*") ^^'^ XN{r) (see the text for the 
definitions). The inset shows the average potential {4>{r)) 
at the Fe atom of the heme (open circles) and its variance 
{{54>{r))'^) (closed diamonds). The dashed lines in the main 
panel connect the points and the solid lines in the inset are 
fits of the simulation data to the function 1 — exp [— (r — a) / A] , 
with A = 12 A and 13 A for the potential and its variance, 
respectively. The calculations have been performed for cytC. 
All results in the main panel are normalized to their corre- 
sponding values at r = 30 A, while the data in the inset are 
normalized to the results obtained for the entire simulation 
box. 



is almost flat, showing virtually no variation of polar- 
ity near the interface, susceptibility Xni^) shows about 
30% increase related to the corresponding increase of the 
interfacial density. This polarity increase is below the 
increment observed at the interface of a non-polar so- 
lute with water )^ reflecting the topological and chemical 
heterogeneity^ii^ of the protein surface. 

Figure [5] suggests that interfacial properties scaling as 
interfacial density will mostly decay to their bulk values 
within ~ 3 hydration layers. However, this rule should 
not be blindly extended to other observables. In partic- 
ular, the electrostatic response of the interface is only 
indirectly affected by density and in principle can have a 
different range of convergence. This is shown in the inset 
in Fig. [5] which presents the accumulation of the average 



and variance of the electrostatic potential 0(r) produced 
by waters from the r-shell at the Fe atom of the heme. 
The accumulation of both the average and the variance 
are much slower than the decay of the interfacial density 
and are in fact comparable in range to the accumulation 
of the cross-correlation Xosir) shown in Fig. U) Specif- 
ically, when these data are fitted to exponential decay 
functions, similarly to what has been done for Xos('")i 
one gets the correlation lengths of 12-13 A. 



B. Dynamical properties 

The long range of the dipolar protein- water correlation 
appearing in susceptibility xos will be masked in xt by 
a larger in magnitude self protein component xoo- The 
two susceptibilities exhibit, however, different dynamics, 
as is shown in Fig. [6] for their loss functions. These are 
calculated from the time correlation function 



(27) 



and the corresponding self and cross correlation functions 

Soa{t) - [((5Mo • 5Ma)]-^ (5Mo(i) • (5M,(0)), (28) 

where a = 0, s. In the case of lys, there is about a fac- 
tor of four difference between the main Debye peaks of 
Xoo(w) (— 14 ns) and Xos(w) (~ 3.5 ns). There is there- 
fore a frequency window in which the two components 
of the overall dipolar susceptibility can be separated. It 
is also worth mentioning here that two relaxation times 
identified here for xos(w), 3.5 ns and 14 ps (lys), bracket 
the typical dielectric J-relaxation band observed between 
protein and water frequency peaks in dielectric loss func- 
tions of protein solutionsi^ The solute-solvent dipolar 
cross-correlations are considered as a plausible cause of 
the (5-dispersion42r— 

Both electrostatic and binding properties of surface 
waters are highly heterogeneous, but this reality is differ- 
ently refiected by the observables. The dynamics of the 
protein dipole moment are essentially single-exponential, 
decaying on the characteristic time of 4-14 ns of protein 
tumbling. The dynamics of the dipole moment of the 
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FIG. 6. Normalized loss function {J^ x" {uj)du} / (jiu]) — 1) 
obtained from the time cross-correlation function So{t) (Eq. 
(|27|l . solid line). The self protein- protein (dash-dotted line, 
p-p) and cross protein-water (dashed line, p-w) loss functions 
are obtained from the corresponding time correlation func- 
tions in Eq. (|28|l . Their relative weights reflect contributions 
of xoo and xoa to Xo ■ The calculations were done for hydrated 
lys. 



first hydration layer are also fairly generic. The corre- 
lation functions decay much faster (see the SM^), but 
still contain a 10-20% slow component with the relax- 
ation time close to that of the protein dipole. This slow 
component can be assigned to waters strongly attached 
to the protein surface4^ 

In contrast to the generic dynamics of the protein and 
first-shell dipoles, the population dynamics of the first 
layer are more specific. The self-correlation functions 
S'Ar(t) of the number of first-shell waters Ni(t) follows 
the rotational dynamics of Mo(i) for cytC and cytB (Fig. 
[7]). It appears that waters in the first layer of these two 
proteins are strongly bound to the protein surface posi- 
tionally, but can rotate relatively freely, resulting in fast 
relaxation of the first-shell dipole moment. On the con- 
trary, waters are weakly bound to the surface of lys and 
ubiq, with significantly faster decays of their populations 
(Fig. El). 

The binding affinity of the first-shell waters is, how- 
ever, heterogeneous. This is seen particularly clear for 
ubiq. Its initial relaxation is two-exponential, with re- 
laxation times of 0.2 ps (58%) and 39 ps (30%). This ini- 
tial fast decay is followed, however, by a plateau indicat- 
ing that about 12% of first-shell waters do not leave the 
protein surface on the simulation time-scale. A similar 
long-time component, contributing to about 17% of the 
correlation function, was found for cytB (see the SM^). 

In addition to the differences in the first-shell popu- 
lation dynamics seen for the two pairs of proteins, the 
first-shell compressibilities {{6Ni)'^) / (Ni) are quite dif- 
ferent for them as well (last column in Table |l|. The 
compressibilities of first shells of lys and ubiq are simi- 
lar to those found for water shells around rigid non-polar 
solutes^i^ while they are much higher for the two cy- 
tochromes. It is not yet clear how this observation re- 
lates to the corresponding differences in the population 



0.8 






CytB 




CytcN 




^ 0.6 








c/?0.4 


Ubiq 






0.2 



~ Lys 











10 



100 

t/ps 



1000 



FIG. 7. Normalized time self-correlation functions of the num- 
ber of waters in the first hydration shell S'jv(t). The data are 
obtained from MD simulations; multi-exponential fits of the 
correlation functions are given in the SM.*^ The dash-dotted 
and dotted lines indicate the self correlation functions Soo(i) 
[Eq. ((25))] of the protein dipole for cytB and cytC, respectively. 



dynamics. 

To connect these observations to our electrostatic prob- 
lem, it seems clear that different observables variably 
report on the protein-water interfacial structure. The 
dipolar response of a protein is particularly sensitive to 
the distribution of the surface charge and might poten- 
tially be considered as a marker of a given protein. High 
variability of eq among the proteins offers a potential for 
applications to protein detection and separation in solu- 
tion. 



V. DIELECTROPHORESIS OF PROTEIN 
SOLUTIONS 

Inhomogeneous external electric field exerts a force on 
the solute dipole i^^i^^ This is given by the expression 



0; 



(29) 



where K = (47r/3)es(xo + Xo"') is the dielectrophoresis 
constant. The solute is dragged toward a stronger field 
if X > (positive dielectrophoresis) or toward a weaker 
field if if < (negative dielectrophoresis) . 

We will use Eq. (|22p to remove dipolar cross- 
correlations and recast K in terms of susceptibilities xoo 
and Xi- One gets 



K 



I 4n 



2, 



1)2/0 



(30) 



where the dipolar density of the solute yo = Ue + 
(47r/3)xoo is introduced in analogy to the dipolar den- 
sity of a polar liquidi^i^ y = {An / 9) /3 pm? ; here, m and 
p are the dipole moment and number density of a liq- 
uid. We have included the component j/e into yo origi- 
nating from electronic polarizability of the solute, which 
has not been considered so far, but needs to be included 
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in high-frequency calculations. This component is ad- 
ditive to the one arising from the permanent dipole^ 
and can be connected to the measurable refractive in- 
dex of the protein no by the Clausius-Mossotti equation, 
Ve = [nl-l)/{nl + 2). 

The two limiting, Lorentz and Maxwell, scenarios are 
worth emphasizing here. In the Lorentz scenario, xi = 
and one gets 



3000 



-ya- 



rn 



The factor in front of yo can be recognized as the Lorentz 
cavity field corrections^ Further, this equation allows 
only positive dielectrophoresis. 

In the Maxwell scenario, xi = Xi^ [Eq. (fT5|) ]. one gets 



A', 



3e. 



1 



M 



2e, 



1 



yo 



2e, 



1 



(32) 



Now, the factor in front of yo is the standard cav- 
ity field of traditional dielectric theories.— Negative di- 
electrophoresis is now allowed when the threshold value 
yo = (^s — l)/(3es) is reached. 

The standard analysis of dielectrophoresis of non-polar 
colloidal suspensiona^i neglects the dipole moment of the 
solute, assuming — in Eq. (j32p . In addition, when a 
dielectric constant eo can be assigned to the material of 
the colloidal particle, in the second summand of Eq. 
(f32|) is replaced by Cs/eo, according to the standard rules 
for the electrostatics of dielectric interfaces ji. The result 
is the traditional constant of dielectrophoresis commonly 
used in the analysis of colloidal suspensions^ 



K = 



£o - £s 
£0 + 2es 



(33) 



This approximation is not applicable to protein solutions 
and we will instead explicitly consider the dipolar re- 
sponse of the protein, with susceptibilities xoo and xi in 
Eq. (pO)) from MD simulations. 

Dielectrophoresis measurements are typically per- 
formed with oscillatory external fields. The static proper- 
ties considered so far need to be replaced with frequency- 
dependent response function according to the standard 
rules of handling time correlation functions.— While the 
frequency-dependent dielectric constant of water is well 
defined and tabulated from laboratory measurements, 
more care is required to define the frequency dependent 
susceptibility Xi('^) a-^d 2/o('^)- 

The standard rules of connecting the response func- 
tions to time correlation functions sugest the following 
form for the frequency-dependent response functions of 
the solute dipole^ 



XOa('^) = XOa 1 - iujSaai~(^) 



(34) 



where a = 0, s and 5oa(<^) is the Fourier-Laplace trans- 
form of the corresponding time correlation functions in 
Eq. (pSj) . These two relations, with Soai^^) obtained from 
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FIG. 8. Re[K(u)], oj = 2nv from Eq. (O with the frequency- 
dependent response functions in Eq. (|34p from MD simula- 
tions and eD(a;) of water from Ref. l54i: the refractive index of 
the protein no = 1.47. 



MD simulations, can be used in Eq. <\22\ to find Xi('^)- In 
doing that, one needs to replace tg with the frequency- 
dependent, complex- valued dielectric constant of water 
€s{uj). Once this is done, one can employ Eq. (|30|) with 
the frequency-dependent susceptibilities and Re[A'(a;)] in 
Eq. (Uni) to calculate the iovceM^ The effect of the elec- 
trolyte can be included into the water dielectric constant 
es(aj) = e£)(aj) -|- <Til{iuj) through the ionic conductivity 
of the solution ai , in addition to the Debye dielectric con- 
stant e£i(aj) = e'^(cj) -|- ie'^{Lo) for the dielectric response 
of pure water. 

The results of calculations of Re[i4r(a;)] for cytB and 
lys are shown in Fig. [8l The dynamics of the protein's 
dipole, Xoo('^) and xos(w) are taken directly from MD 
simulations (see the SM^ for the fits of the relaxation 
functions and the list of the relaxation times). The 
frequency-dependent dielectric constant of water eoi^) 
is taken from Ref. [13 and no = 1.47 is assigned to 
the protein. The electrolyte conductivity in the range 
Ui ~ 1 — 100 mSm^^ typically employed in the labora- 
tory measurements^! does not affect the outcome of the 
calculations. 

The main result of these calculations is a dominance of 
yo(uj), arising from the protein permanent dipole, in the 
overall dielectrophoresis response. This result holds for 
all charged proteins studied here. In contrast, even static 
dielectrophoresis constant is negative for ubiq, suggesting 
negative dielectrophoresis in the entire frequency range. 
Numerical difficulties in obtaining Sos{t) for ubiq have 
prevented us from presenting K{lu) for this protein. 



VI. THZ ABSORPTION OF PROTEIN 
SOLUTIONS 

Positive xij instead of the negative values in both the 
Lorentz and Maxwell scenarios, alter the theory predic- 
tions concerning the absorption of THz radiation by pro- 
tein solutions. We have recently suggested a model^ 
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based on the formalism of response functions discussed 
in Sec. 1111) to address this problem. The study was mo- 
tivated by experimental reports of the positive slope of 
the absorbance increment with an increasing concentra- 
tion of the protein in solution, following by a non-linear 
downward turn of the concentration dependence^>^ (in- 
set in Fig. O. Such a trend, recorded in the 1-3 THz 
frequency range, contradicts the traditional view that 
adding a protein, less polar than water, should lower 
the solution polarity and the corresponding radiation ab- 
sorbance. Indeed, the model calculations in Ref. H could 
not account for the observations without demanding a 
significant increase of the effective dipole of the protein- 
water interface. That development was, however, based 
on the Maxwell scenario for the interface susceptibility 
and thus negative interface dipole. Since the present for- 
malism goes beyond the Maxwell picture and the numer- 
ical results allow a positive interface dipole, this problem 
needs to be revisited. 

We are looking at the absorbance of the electromag- 
netic radiation^ 



(35) 



where c is the speed of light and the superscript "T" em- 
phasizes that we are considering the susceptibility of the 
solution x(<^) in the direction of the electric field perpen- 
dicular (transversal) to the direction of light propagation 
determined by the wave- vector. The property reported 
experimentally is the relative increment Aa(w)/as(cj) = 
Q!mix(w)/as(cij) — 1 of the solution (mixture) absorbance 
over the absorbance of the homogeneous liquid as(w). 

The derivation of the solution absorbance directly fol- 
lows from the generalized response function in Eq. pip . 
The derivation steps follow Ref. and are briefly sum- 
marized in Appendix [B1 One gets the following result for 
the relative susceptibility increment 



Ax' {^)/xi (^) = -^0 

"e,,(w) + 2 Stt 



47r 

1 - Yes(w)xi(w) 



2/o(w)??o 



where Ax'^(a;) 



es(iLj) - 1 



Xi{u)es{uj)Io{rio,R) 



(36) 



Xmix('^)-xnw) and x^A^) = (e,s(w)- 



l)/(47r). Further, 770 is the volume fraction of the solutes 
and /q {"Hq I R) represents the mutual polarization of the 
solutes by their permanent dipoles aligned by the field of 
radiation, a non-ideal solution effect. It is given by Eq. 
(|B7p in Appendix [Bl and has been tabulated as a function 
of ?7o and R in Ref. y assuming hard-sphere structure fac- 
tor for the solutes in solution. This latter approximation 
has been used in the calculations presented in Fig. [51 

All the parameters in Eq. ([36)) are defined in our calcu- 
lations of the dielectrophoresis response and are directly 
applied to the calculations shown in Fig. [51 for cytB. The 
solid line shows lS.a{ijj) / as{uj) calculated at 2.25 THz of 
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FIG. 9. Absorbance increment calculated at the frequency 
of 2.25 THz vs the molar fraction r^o of cytB in water. The 
calculation with Xi('^) from MD simulations (solid line) is 
compared to Xi('^) from either Maxwell or Lorentz scenarios 
(dashed line, labeled "M-l-L") which give very close results 
in this frequency range. The dash-dotted line is the calcu- 
lation performed with the restriction (in /3)es{ijj)xi{^) ~ 1 
corresponding to zero xos in Eq. (|22p . The bold dashed line 
is arbitrary drawn to illustrate a crossover between two Un- 
ear slopes when proteins in solution start to affect each other 
and prevent full convergence of xos {r) characterized by a long 
correlation length of ~ 2.4 nm [FigJ31. The inset shows ex- 
perimental data reported in Ref. 53 for a five-helix bundle 
protein Ag_g5 at two temperatures indicated in the plot. The 
dotted lines in the inset connect the points. 



radiation used in the experiment)^ The result is an ob- 
viously positive slope of the absorption increment, which 
arises from both the interface and permanent dipoles 
re-enforcing each other. In contrast, both Lorentz and 
Maxwell scenarios suggest a negative slope (dashed line 
in Fig. [9]). 

The non-linear dependence on the solute volume frac- 
tion present in /o(i?, ^^o) in Eq. ()36p turns out to be insuf- 
ficient to bend the concentration dependence downward 
(inset in Fig. [9]). Several possible long-range effects are 
still missing from our analysis. The structure factor of 
solvated proteins is estimated from its hard-sphere limit 
and can be modified by long-ranged interactions. In addi- 
tion, correlation between proteins' dipole moments, rep- 
resented by the Kirkwood factor, has not been included 
in the present calculationi^ We, however, want discuss 
yet another possibility related to the long decay of the 
solute-solvent dipolar correlations shown in Fig. [H 

The correlation length of A ~ 2.4 nm found in Fig. [4l 
suggests that susceptibility xos becomes affected when 
proteins in solution come closer than 2A, which roughly 
corresponds to the position of the peak in the experi- 
mental concentration dependence in the inset of Fig. [9l 
If mutual effect of the proteins in solution docs not allow 
Xos to saturate, it implies that, according to Eq. ([22|l . 
Xi should change with the concentration and approach 
the limit (47r/3)es(cj)xi(i^) — 1. This limit for the ab- 
sorbance is shown by the dash-dotted line in Fig. [51 The 
overall concentration dependence is expected to exhibit a 
crossover from one linear slope to another in the concen- 
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tration range where proteins start to afFect each other, as 
is schematically shown by the bold dashed line in Fig. [9l 
The experiment shows a qualitatively similar crossover 
behavior. Further, no cross-over of the concentration 
dependence of the absorption coefficient was found for 
solutions of aminoacids.— This observation suggests a 
much shorter length-scale of dipolar cross-correlations for 
hydratcd aminoacids compared to proteins. The reason 
might be related to a more homogeneous structure of wa- 
ters around aminoacids and thus a higher contribution of 
the first hydration layer to xqsj which is almost absent 
for proteins (Fig. 

A brief comment on the experimental temperature de- 
pendence of the absorption shown in the inset of Fig. [S] is 
relevant here. The permanent dipole susceptibility xoo is 
proportional to /?. The slope of the absorbance with con- 
centration is therefore expected to increase with lowering 
temperature, as is indeed qualitatively observed. 

These calculations and qualitative comparisons with 
experiment suggest that most of the polar response of 
the protein charges is frozen on the time-scale of THz 
radiation, which therefore allows one to probe the po- 
larization of the protein- water interface projected on the 
dipole Afg"*. This interface polarization is distinctly dif- 
ferent from the Lorentz-Maxwell scenario and this fact 
is reflected in the positive slope of the absorbance with 
increasing concentration. 



VII. SUMMARY 

In conclusion, we have utilized long simulation trajec- 
tories and large system sizes to systematically study the 
dipolar susceptibility [Eq. ([7])] of hydrated proteins. This 
property exhibits large variation among the proteins, sug- 
gesting its possible use as a "marker" of a protein in so- 
lution. The large variance in the protein susceptibility 
also opens the door to their detection and separation in 
solution, in particular by dielectrophoresis of the protein 
solutions. 

The permanent dipole susceptibility is a combination 
of the variance of the intrinsic protein dipole and a cross- 
correlation between the protein and water dipoles. The 
cross-correlation component is found to be negative for 
uncharged ubiq, but is positive for the charged proteins. 
It is also long-ranged, decaying into the bulk on the cor- 
relation length of about 2 nm. A similar, but somewhat 
smaller, decay length is found for the electrostatic poten- 
tial produced by the hydration water inside the protein. 

Dipolar solute-water cross-correlations are significant 
and cannot be neglected. In the case of ubiq the self and 
cross correlation components nearly cancel each other, 
but produce a negative overall susceptibility. This pro- 
tein therefore demonstrates negative dielectrophoresis 
and, correspondingly, a dia-electric dipolar response. 

Thermodynamic stability prohibits dia-electric re- 
sponse (eg < 1) for bulk dielectricsii This limitation 
docs not apply to a finite-length polar response;^ and in 



fact the wave- vector dependent dielectric constant es(fc) 
is negative in a certain range of fc-values for most po- 
lar liquidsi^ Thermodynamic arguments also do not rule 
out a dia-electric response of solutes in solution. Nega- 
tive dielectrophoresis is in fact fairly common for colloidal 
suspensions.— Nevertheless, this is the first report of neg- 
ative dielectrophoresis of proteins by numerical simula- 
tions. 

An important question posed by the present study and 
requiring further investigation is whether negative values 
of xo found here for ubiq are general for neutral pro- 
teins. If the answer to this question is affirmative, dielec- 
trophoresis of hydrated proteins is expected to be sensi- 
tive to the buffer pH (as indeed has been reported^l). In 
addition, negative dielectrophoresis, and thus dia-electric 
response, should be common when the buffer pH ap- 
proaches protein's isoelectric point. Future simulations 
and laboratory measurements are required to shed more 
light on this intriguing perspective. In terms of phys- 
iological conditions of protein activity, this possibility 
would imply a range of pH values in which a protein is 
nearly insensitive to significant gradients of local electric 
fields, i.e., "invisible" to local fields. 

The interface polarization is a small portion of the 
static or low-frequency dipolar response, but gains in im- 
portance at high frequencies, above the protein Debye 
peak representing its rotation in solution. In this fre- 
quency range, the polar response of the protein charges 
is dynamically arrested and the response of much faster 
interfacial waters starts to show up. Here, different di- 
electric scenarios (such as Lorentz or Maxwell recipes) 
predict distinctly different outputs. Numerical simula- 
tions presented here show a broad range of possibilities, 
beyond those emphasized by these two electrostatic lim- 
its, depending on the distribution of the protein sur- 
face charge. These different outcomes are reflected by 
the THz absorbance of the protein solution. In particu- 
lar, the sign of the slope of the THz absorbance vs the 
protein concentration reports on whether the interface 
dipole is oriented along the polarizing held or opposite 
to it. A parallel orientation corresponds to a positive 
slope, whereas an anti-parallel orientation (such as in 
the Maxwell scenario) yields a negative slope. The mea- 
surements of the THz absorbance give therefore a direct 
access to the polarization pattern of the protein-water 
interface. 
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Appendix A: Derivation of Eq. (|22p 

The function x(ki, k2), first calculated in Refs.[36l and 
l60l . describes the response of a polar liquid outside a 
void. The latter is defined by the step function 0o(k) 
[Eq. (PT|) ]. This function and its conjugate 9 = 1 — ^o, 
l(ki,k2) = (27r)'^(5(ki — k2) together form a set of or- 
thogonal functions projecting the polar response and cor- 
responding electric fields inside and outside the solute. 
They obey the following orthogonality and multiplica- 
tion rules: 6 * 6q = 0, 6 * 6 ~ 9, Oq * Oq ~ Oq, where, as 
above, asterisks refers to the k-space convolution. 

The inhomogeneous dipolar response function can then 
be written as follows 



Xs - axs *0o*G ^ *0o*Xs 



(Al) 



Here, Xs(ki,k2) = Xs(ki)l(ki, k2) and G = 9o*Xs*do, 
which in terms of explicit wave- vectors implies 

G(ki,k2) = ^o(ki -k') *x.(k',k") *0~o(k",k2). (A2) 

Equation (jAl[) with a = 1 presents an exact solution 
of the problem of dipolar response of a polar liquid inter- 
facing a spherical void j'^^'^" It does not, however, antici- 
pate the formation of the specific orientational structure 
at the interface characteristic of liquid water 4^^^!^ There- 
fore, a = Xi/xY 7^ 1 is introduced in Eq. (|Aip to account 
for possible deviations from the Maxwell scenario of the 
surface polarization. 

From the definition of the response function, the pro- 
jection of the polarization inside the solute is 



9o * P = (1 - a)0o * Xs * Eo. 



(A3) 



The polarization projection is of course zero in case 
of a = 1 since this is how the response function was 
constructed.— A modification of the inhomogeneous re- 
sponse by a 7^ 1 in the second term in Eq. (jAl[) is 
physically equivalent to creating a non-zero dipole as- 
sociated with the solute. In case of a uniform external 
field Eo(k) = So ^^zEq this dipole becomes 



Oz 



z • 0o(-k) * P(k) = (1 - a)x^iO)noEo. (A4) 



When the electric field is produced by the solute dipole 
[Eq. (dH)], this field can be written as the projection 
of the dipolar field outside the solute excluded volume 
Eq = 0*Eci. Here, E^; is the Fourier-space electric field of 
a point dipole, without the space cutoff equating the field 
to zero inside the solute and responsible for the appear- 
ance of the spherical Bessel function in Eq. (fTO)) . With 
this representation one gets from Eq. (jA3p 



(A5) 



9o * P = (1 - a)eo * X. * * Ed 



If XsC^) here is replaced with its continuum limit Xs(0), 
the integral becomes identically zero because of the or- 
thogonality of ^0 and 9. The result is = P(0) in Eq. 

(Ha- 



lf we now choose the direction of z-axis along Mo and 
take the projection of the solute dipole on z, we get 



M,,,/A/o = z-x(0,k)*T(k).z. 



(A6) 



A note on how to correctly take the k ^ limit in the 
k-space tensor functions is appropriate here. Since both 
z and k impose axial symmetry on the liquid, which is 
isotropic in the direct space, one has to take k parallel 
to z when taking the fc — > limit to avoid imposing a bi- 
axial symmetry. In this prescription, Tzz{0) = — (Stt/S) 
and one gets for the homogeneous component of x(0,k) 
(first summand in Eq. (|Aip ) 



z-x.(0)-T(0)-z 



2(e.-l) 
3ec 



(A7) 



Similarly, one can write down the inhomogeneous com- 
ponent of the polar response from Eq. (fTTj) as 



(xfr'z-Xo(0,k)*T(k).z = 



ttR^xHO) 



dfcii(fci?)2(z.x,,(k) - D-z), 



(A8) 



where {. . refers to the average over the solid angles 
of the unit vector k and 

JiikR) 



9oik)/no 



kR 



(A9) 

and (z • • z)(jj. = (2/3) one gets in the continuum limit 
A: — !■ for the response function Xs(k) 



has been used. Using the relations (z • • z)^ 



(xf)"'2-Xo(0,k)*T(k).z 
where 



-(87r/9)(e, - 1), (AlO) 



{6R/tt) ji{kRfdk^l (All) 

and Eqs. (fT2|) for the A: = values of Xs '"^(0) have been 
used. Combining Eqs. (jA7P and (jAlOp . one arrives at Eq. 

A similar procedure can be applied to calculate the 
dipole moment of the solvent in uniform external field 
Eo(k) ~ (5o,kzi?o- The dipole moment in this case is 



z.0~(-k')*x(k',k)*Eo(k). 



(A12) 



According to the preceding arguments, ^*X = X— (1 — 
* Xs- From this relation one gets 

Ms,. - Miiq + = A/iiq - Mo, + (A13) 

where Mq. is given by Eq. (jA4p and Afuq is the dipole in- 
duced by Eq in a homogeneous liquid without the solute. 
Further, Mq^ in Eq. (|A13|) is the "cavity dipole" associ- 
ated with inserting the excluded solute volume into the 
polarized liquid. It is given as 

Mo=J(^2oi?o) = XI - xf'(O) (1 - Xi/xf) . (A14) 

Equations (|A13|) and (|A14p yield Eq. Q for the overall 
solution dipole. 
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Appendix B: Derivation of Eq. (|36|) 

We want to calculate the dipole moment M"^{u:) of 
the solution induced by a polarized electromagnetic wave 
propagating along the z-axis of the laboratory frame and 
having the electric field vector along the a;-axis of the 
same frame. The wave- vector is therefore along z and 
the response is transversal. The susceptibility in Eq. (|35|) 
becomes 



where 



(Bl) 



The formalism of response functions^i^ described in 
Appendix [X] gives the following prescription for the cal- 
culation of the dipole moment of the solution projected 
on X 



Xoo^oNoEo + X • ^(k) * x(k, k') * Eo(k'). (B2) 



The first term in this equation (dependencies on fre- 
quency are omitted for brevity) is the direct polarization 
of the permanent dipoles of Nq solutes in the mixture 
and the second term is the polarization of the solvent. 

The formalism of Appendix |X] is now extended to an 
ensemble of solutes. This extension is straightforward if 
polarization fields of the liquid around different solutes 
are uncorrelated.~ From Eqs. (fTT|) and (jAip one gets 



X(ki,k2) = Xs(ki)(5ki,k2 -aX]^"(''i'^2)e* 



k2 -r^ 



(B3) 



where the sum runs over the Nq solutes in the solution. 
In addition, the source of the electric field now includes 
the homogeneous external field and the field of all dipolar 
solutes polarized by it 



Eo(k) = -Bo'^o.kx + Xoa^oEa J2 T(k) ■ xe 



(B4) 



According to the rules of calculating the k-space con- 
volutions explained in Appendix \^ the substitution of 
the uniform external field (first summand in Eq. (jB4p ) 
into Eq. (|B2|) leads to 



VEo 



1 - 770 1 



2e., 



(B5) 



where ?7o = ^o^o/^ is the solute volume fraction. 

Similarly, the use of the second, dipolar summand from 
Eq. dm in Eq. (|B2| yields 



^2 47r J, 

= Y'7oXs (O)Xoo 



1 - aIoim,R) 



2(e. - 1) 



263 + 1 

(B6) 

Here, one employs Txx{0) = (47r)/3 for the transverse 
component of the dipolar tensor and 

QR /"CO 

loim.R)^— 3i{kR?SQ{k)dk, (B7) 
I" Jo 



(B8) 



is the density structure factor— of the solutes in the so- 
lution. 

The integral /q (770 , R) is equal to unity for an ideal 
solution with S'o(fc) = 1 [as in the case of Eq. (jAlip j. 
It is also tabulated in Ref. in terms of a polynomial 
interpolation in the solute size and volume fraction in 
the hard-sphere approximation for So{k). This latter ap- 
proximation was used in the calculations shown in Fig. [9l 
Finally, Af ^ = xoo^oNqEq + Mf + from Eqs. ((B5|) 
and (jB6[) yields the susceptibility increment in Eq. ([36]). 

We note that ye from protein polarizability has been 
added to (47r/3)xoo in Eq. ([5^1) . One can use experimen- 
tal refractive index of the protein powder in the same fre- 
quency range as for the solution absorption to estimate 
the contribution of the protein intrinsic dipole from the 
Clausius-Mossotti equation {no{uj)^ — l)/{no{uj)^ + 2) = 
ye + (47r/3)xoo the use of which should be restricted to 
THz and higher frequencies. 
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